function heat_capacity_versus_collection()
    n = 40;
    p = 0;
    b = 0;
    isingmat0 = isingmat_ini(n, p);
    flip_num = 1;
    T = 2.25;
    
    step_max = 5000;
    num = 100;
    collect_step = linspace(1, (step_max - step_max / num + 1), num);
    
    steps = 20000;
    sam_len = 1500;
    
    Clist = zeros(1,num);
    
    for cy = 1:steps
        [trans_mat, km, mm]= trans_state(isingmat0, n, T, b, flip_num);
        isingmat0 = trans_mat;
    end
    
    for i = 1:num
        step = collect_step(i);
        en_list = zeros(1,sam_len);
        for j = 1:sam_len
            isingmat = isingmat0;
            for k = 1:step
                [trans_mat, mk, mm]= trans_state(isingmat, n, T, b, flip_num);
                isingmat = trans_mat;
            end
            en = ising_energy_in(isingmat, n, b);
            en_list(j) = en;
        end
        
        C = var(en_list)/(T*T); 
        Clist(i) = C;

    end

    figure(1)
    scatter(collect_step, Clist)
    xlabel('collect-step');
    ylabel('C (a.u.)');
    title_string = sprintf ('C versus collect-step (n = %d, T = %.1f, b = %.1f, p = %.1f)', n, T, b, p);
    title(title_string);
    
    
end


function E = ising_energy_in(mat, n, b)
    J = 1;
    E = - b * sum(sum(mat));
    for i = 1 : n
        for j = 1 : n
            if(i == 1)
                if(j==1)
                    sumi = mat(1, 2) + mat(2, 1) + mat(1, n) + mat(n, 1);
                elseif(j==n)
                    sumi = mat(1, n-1) + mat(2, n) + mat(n, n) + mat(1, 1);
                else
                    sumi = mat(1, j-1) + mat(1, j+1) + mat(2, j) + mat(n, j);
                end
            elseif(i == n)
                if(j==1)
                    sumi = mat(n, 2) + mat(n-1, 1) + mat(1, 1) + mat(n, n);
                elseif(j==n)
                    sumi = mat(n, n-1) + mat(n-1, n) + mat(1, n) + mat(n, 1);
                else
                    sumi = mat(n, j-1) + mat(n, j+1) + mat(n-1, j) + mat(1, j);
                end
            else
                if(j==1)
                    sumi = mat(i-1, 1) + mat(i+1, 1) + mat(i, 2) + mat(i, n);
                elseif(j==n)
                    sumi = mat(i-1, n) + mat(i+1, n) + mat(i, n-1) + mat(i, 1);
                else
                    sumi = mat(i-1, j) + mat(i+1, j) + mat(i, j-1) + mat(i, j+1);
                end
            end
            E = E - 0.5 * J * mat(i,j) * sumi;
        end
    end
    return
end
